package heli;

//Author Hannes Hassler
//


public class SunCalculations {
	
	/**
	 * 
	 * Author Hannes Hassler, (C) 2011
	 * license: Eclipse Public License - v 1.0
	 *          (http://www.eclipse.org/legal/epl-v10.html)
	 * 
	 * Approximation Formulas to calculate the Azimuth and elevation of the sun
	 * according to 
	 * http://de.wikipedia.org/wiki/Sonnenstand	  
	 * http://en.wikipedia.org/wiki/Julian_day
	 *              
	 * For better verification I maintained in the method names also
	 * the short hand notation from the Wikipage, combined with the actual meaning.
	 * The Wiki-Author(s) point out that these calculations follow 
	 * certain simplifications, citing more exact calculations, too. (VSOP87)
	 * The Wiki-Author(s) claim these approximations precise whithin an accuracy of 0.01�
	 * and years 1950 to 2050.
	 *  
	 * 
	*/
	    
	  /** Julian date, Ref JD*/
	  public static double toJulian_JD(int year, int month, int day, double hour, double minute, double second) {
		  //Taken from http://en.wikipedia.org/wiki/Julian_day; Calculation
		  double jd;		  
		  int a = (14-month)/12;
		  int y = year +4800-a;
		  int m = month +12*a-3;		  
		  int jdn = day + (153*m+2)/5+365*y+ y/4 -y/100+y/400-32045;		  
		  jd = jdn + (hour-12)/24 + (minute/1440)+ second/86400;		  
		  return jd;
	  }
	  
	  /**days since 1.Jan.2000, Ref n*/
	  public static double days_since_J2000_n(double jd) {
		  return jd-2451545.0;
	  }

	  /**medium ecliptian length (raw); Ref L_raw*/
	  public static double med_ecl_L_raw(double n) {
		  return 280.46+0.9856474*n;
	  }
	  
	 /**medium ecliptian length; Ref L (multiples of 360 removed from L_raw)*/
	  public static double med_ecl_L(double n) {
		  double rawVal=med_ecl_L_raw(n);		  
		  return remove_multiples(rawVal,360);
	  }
	  
	  /**medium anomaly(raw); Ref g_raw*/
	  public static double med_anomaly_g_raw(double n) {
		  return 357.528+0.9856003*n;
	  }
	  
	  private static double remove_multiples(double val2Remove, int modulus) {
		  int overShoots=0;
		  if (val2Remove>=0) {
			  overShoots=(int) val2Remove/modulus;
		  } else {
			  overShoots=((int)-val2Remove/modulus)-1;
		  }
		  return val2Remove-overShoots*modulus;
	  }
	  
	  private static double toPositiveQuadrant(double val2Positive) {		  
		  double rawVal = remove_multiples(val2Positive,360);		  
		  if (rawVal>180) rawVal-=180;		  
		  return rawVal;
	  }
	  
	  /**medium anomaly(raw); Ref g*/
	  public static double med_anomaly_g(double n) {		  
		  return remove_multiples(med_anomaly_g_raw(n),360); 
	 }

	  /**medium ecliptic length; Ref Lambda */
	  public static double med_ecl_length_Lambda(double g, double L) {
		  //assuming numerical eccentricity e=0.0167 c.Wikipage
		  double g_rad=toRadians(g);
		  double rawVal=L+1.915*sin(g_rad)+0.02*sin(2*g_rad);
		  return remove_multiples(rawVal,360); 
	  }
	  
	  /** ecliptic; Ref epsilon*/
	  public static double ecliptic_epsilon(double n) {
		  return remove_multiples(23.439-0.0000004*n,360);
	  }
	  
	  /** right ascension; Ref alpha*/
	  public static double right_ascension(double epsilon, double Lambda) {
		  double epsilon_r=toRadians(epsilon);
		  double Lambda_r=toRadians(Lambda);
		  double rawVal = toDegrees(atan(cos(epsilon_r)*tan(Lambda_r)));
		  return toPositiveQuadrant(rawVal);
	  }
	  
	  /**declination; Ref delta*/
	  public static double declination_delta(double epsilon, double Lambda) {
		  double epsilonRad=toRadians(epsilon);
		  double LambdaRad=toRadians(Lambda);
		  double deltaRad=asin(sin(epsilonRad)*sin(LambdaRad));
		  return toDegrees(deltaRad);
	  }
	  
	  /**T0*/
	  public static double T0(double jd0) {	 
		  return (jd0-2451545.0)/36525; 
	  }
	  
	  /**average sidereal time (hourly, raw); Ref theta__h_G*/
	  public static double average_sidereal_time_hourly_theta__h_G_raw(double T0, int hour, int minute, int second) {
		  double T = hour+((double)minute)/60+((double) second)/3600;
		  return 6.697376 + 2400.05134*T0 + 1.002738*T;
	  }
	  
	  /**average sidereal time (hourly); Ref theta__h_G*/
	  public static double average_sidereal_time_hourly_theta__h_G(double T0, int hour, int minute, int second) {
		  double rawVal=average_sidereal_time_hourly_theta__h_G_raw(T0, hour, minute, second);
		  return remove_multiples(rawVal,24); 
	  }
	  
	  /**hour angle of equinox; Ref theta*/
	  public static double hour_angle_of_equinox_theta(double theta__h_G, double lambda) {
		  double theta_G=theta__h_G*15;
		  return theta_G+lambda;		  
	  }
	  
	  /**hour angle of sun for location; Ref tau*/
	  public static double hour_angle_of_location_tau(double theta, double alpha) {
		  return theta-alpha;
	  }
	  
	  /**Azimut; Ref A*/
	  public static double azimuth_A(double tau, double phi, double delta) {
		  double tauRad=toRadians(tau);
		  double phiRad=toRadians(phi);	
		  double deltaRad=toRadians(delta);		  
		  double denominator = (cos(tauRad)*sin(phiRad))-
				               (tan(deltaRad)*cos(phiRad));	        	
		  double azimuthRad=atan(sin(tauRad)/denominator);
		  double azimuthDeg=toDegrees(azimuthRad);
		  if (denominator<0) azimuthDeg+=180;
		  if (azimuthDeg>180) azimuthDeg-=360;
		  return azimuthDeg;
	  }
	  
	  public static double azimuth(
			    //date, time
			    int year,
				int month,
				int day,
				int hour,
				int minute,
				int second,
				
				//geocodes 
				double phi, //N
				double lambda //O
			  ) {
		  
		    double jd=SunCalculations.toJulian_JD(year, month, day, hour, minute, second);
			double jd0=SunCalculations.toJulian_JD(year, month, day, 0, 0, 0);	
			double n = SunCalculations.days_since_J2000_n(jd);	
			double g = SunCalculations.med_anomaly_g(n);
			double L = SunCalculations.med_ecl_L(n);
			double Lambda = SunCalculations.med_ecl_length_Lambda(g, L);
			double epsilon = SunCalculations.ecliptic_epsilon(n);
			double alpha = SunCalculations.right_ascension(epsilon, Lambda);
			double delta = SunCalculations.declination_delta(epsilon, Lambda);
			double T0 = SunCalculations.T0(jd0);
			double theta__h_G = SunCalculations.average_sidereal_time_hourly_theta__h_G(T0, hour, minute, second);
			double theta = SunCalculations.hour_angle_of_equinox_theta(theta__h_G, lambda);
			double tau = SunCalculations.hour_angle_of_location_tau(theta, alpha);
			double azimuth_A = SunCalculations.azimuth_A(tau, phi, delta);
			
		  return azimuth_A;
		  
	  }
	  
	  public static double elevation(
			    //date, time
			    int year,
				int month,
				int day,
				int hour,
				int minute,
				int second,		
				
				//geocodes 
				double phi, //N
				double lambda //O
				) {
		  
		    double jd=SunCalculations.toJulian_JD(year, month, day, hour, minute, second);
			double jd0=SunCalculations.toJulian_JD(year, month, day, 0, 0, 0);	
			double n = SunCalculations.days_since_J2000_n(jd);	
			double g = SunCalculations.med_anomaly_g(n);
			double L = SunCalculations.med_ecl_L(n);
			double Lambda = SunCalculations.med_ecl_length_Lambda(g, L);
			double epsilon = SunCalculations.ecliptic_epsilon(n);
			double alpha = SunCalculations.right_ascension(epsilon, Lambda);
			double delta = SunCalculations.declination_delta(epsilon, Lambda);
			double T0 = SunCalculations.T0(jd0);
			double theta__h_G = SunCalculations.average_sidereal_time_hourly_theta__h_G(T0, hour, minute, second);
			double theta = SunCalculations.hour_angle_of_equinox_theta(theta__h_G, lambda);
			double tau = SunCalculations.hour_angle_of_location_tau(theta, alpha);
			
			double elevation_h = SunCalculations.elevation_h(tau, phi, delta);
			double elevation_h_R = SunCalculations.elevation_refrCorrection_h_R(elevation_h);
		  
		  return elevation_h_R;
      }
	  
	  /**Elevation; Ref h*/
	  public static double elevation_h(double tau, double phi, double delta) {
		  double tauRad=toRadians(tau);
		  double phiRad=toRadians(phi);	
		  double deltaRad=toRadians(delta);	
		  double elevationRad=asin(
				               cos(deltaRad)*cos(tauRad)*cos(phiRad)+
				               sin(deltaRad)*sin(phiRad)
				               );		  
		  return toDegrees(elevationRad);
	  }
	  
	  /**Elevation refraction correction; Ref h_R*/
	  public static double elevation_refrCorrection_h_R(double h) {		  
		  //Assuming atmospheric pressure 1010 millibar, 
		  //and Temperature 10°C
		  double denominatorDeg=h+(10.3/(h+5.11));
		  double denominatorRad=toRadians(denominatorDeg);
		  double R=1.02/tan(denominatorRad); 
		  double h_R=h+R/60;		               		  
		  return h_R;
	  }
	  
	  
	  private static double toRadians(double angDeg) {
		  return java.lang.Math.toRadians(angDeg);
	  }
	  
	  private static double toDegrees(double angRad) {
		  return java.lang.Math.toDegrees(angRad);
	  }
	  
	  private static double tan(double ang) {
		  return java.lang.Math.tan(ang);
	  }
	  
	  private static double atan(double ang) {
		  return java.lang.Math.atan(ang);
	  }
	  
	  private static double cos(double ang) {
		  return java.lang.Math.cos(ang);
	  }
	  
	  private static double sin(double ang) {
		  return java.lang.Math.sin(ang);
	  }
	  
	  private static double asin(double ang) {
		  return java.lang.Math.asin(ang);
	  }
	  
}
